arXiv:1502.03662vl [q-bio.NC] 12 Feb 2015 


A Rich Spectrum of Neural Field Dynamics in the Presence of 
Short-Term Synaptic Depression 

He Wang,^ Kin Lam,^ C. C. Alan Fung,^ K. Y. Michael Wong,^ and Si Wu^ 

^Department of Physies, Hong Kong University of Seienee and Teehnology, Hong Kong, China 
State Key Laboratory of Cognitive Neuroseienee and Learning, 

IDC/MeCovern Institute for Brain Researeh, Beijing Normal University, Beijing 100875, China 

(Dated: February 13, 2015) 

In continuous attractor neural networks (CANNs), spatially continuous information such as ori¬ 
entation, head direction, and spatial location is represented by Gaussian-like tuning curves that can 
be displaced continuously in the space of the preferred stimuli of the neurons. We investigate how 
short-term synaptic depression (STD) can reshape the intrinsic dynamics of the GANN model and 
its responses to a single static input. In particular, GANNs with STD can support various complex 
firing patterns and chaotic behaviors. These chaotic behaviors have the potential to encode various 
stimuli in the neuronal system. 
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I. INTRODUCTION 

To understand the brain’s operation, it is important to 
consider the range of firing patterns and their conditions 
of occurrence, which will shed light on how information 
is encoded in the brain [1]. In the processing of contin¬ 
uous information such as object orientation and spatial 
location, firing patterns are found to be localized in the 
space of preferred stimuli of the neurons, normally tak¬ 
ing up a Gaussian-like profile mm- Thus an interesting 
question is whether these profiles are stable in time and 
in space and, if not, what other dynamical states will 
replace them. 

Gaussian-like profiles play an important role in both 
experiments and theory. The tuning curves, i.e., the func¬ 
tional dependence of the neuronal response on the inputs 
and the preferred stimuli of the neurons, are observed as 
Gaussian-like profiles in various animal experiments. For 
example, head direction cells in anterior thalamus and 
postsubiculum were found in rodents Hi]. Their tuning 
curves are Gaussian-like and centered at their preferred 
head directions. Place cells in hippocampus are another 
example HIHI. Place cells’ activities are observed to de¬ 
pend on the location of the animal within the environ¬ 
ment. The tuning curve is also a Gaussian-like function. 
Moving direction cells in middle temporal (MT or V5) 
cortex in macaques are found to be selective to object 
moving directions, which also have Gaussian-like tuning 
curves [9l[TO]. 

In neural field models processing continuous informa¬ 
tion, Gaussian-like tuning curves are steady states of the 
network dynamics, and remain stable when their posi¬ 
tions are displaced in the space of the preferred stimuli 
of the neurons. These neural field models are called con¬ 
tinuous attractor neural networks (CANNs), since the 
Gaussian-like tuning curves are attractors of the network 
dynamics. Recent evidence supporting the existence of 
continuous attractors was reported by Wimmer et al, 
who discovered activity anticorrelations on the opposite 


sides of tuning curves as predicted by the CANN model 
in the prefrontal cortex of monkeys [nHi3]. 

The ability to support these bump attractors is effected 
by couplings between neurons in CANNs. However, in 
reality, couplings between neurons are not quenched. 
They depend on firing histories of presynaptic neurons. 
Tsodyks et al. found that synaptic efficacy decreases 
with firing history [niiiii. Furthermore, they proposed 
that this decline in synaptic efficacy is due to the slow 
dynamics of the recovery process of neurotransmitters. 
The recovery of neurotransmitters is of the order of 100 
ms. This short-term decline in synaptic efficacy is called 
short-term synaptic depression (STD). 

Various effects of STD on CANN have been studied 
in the literature. For instance, Fung et al reported that 
CANN with STD can support four different phases ac¬ 
cording to strengths of inhibition and STD [16] . In par¬ 
ticular, STD can drive traveling bumps in the attractor 
space (see Fig. [^c)). These moving profiles happen in 
the absence of external inputs. Additionally, York and 
van Rossum reported a similar result, but with a uniform 
background current m- 

Moreover, the network’s response can be modulated by 
the interplay between the STD-driven intrinsic dynamics 
and a moving stimulus. Without STD, a bump-shaped 
neuronal activity profile tracks a moving stimulus with 
a delay. However, with a proper strength of STD, the 
network activity profile can move ahead of the moving 
stimulus. Effectively, the network activity profile is lo¬ 
cated at a future position of the moving stimulus. It 
can be used to implement an anticipation mechanism, 
which can compensate inherent delays in the neural sys¬ 
tem, thus achieving real-time tracking m- 

For a static input, short-term synaptic depression 
can also drive periodic excitements of neuronal activity 
nsj. We proposed that these periodic excitements enable 
CANNs to support representations of multiple stimuli al¬ 
most simultaneously [20]. These periodic excitements 
provide a plausible mechanism for resolution enhance¬ 
ment in transparent motion uniEi]. 
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FIG. 1. (Color online) Intrinsic Behaviors. The color scale 
shows U{x,t). (a) Silent state. Parameters: k = 0.8, /3 = 0.2. 
(b) Static bump. Parameters: k = 0.8, 13 = 0.005. (c) Moving 
bump. Parameters: k = 0.8, [3 = 0.05. (d) Uniform firing. 
Parameters: /c = 1 x 10“^, [3 = 0.02. (e) Homogeneous spikes. 
Parameters: k = 1 x 10“^, [3 = 0.023. (f) Spikes and anti¬ 
spikes. Parameters: k = lx 10“^, /3 = 0.0245. For all (a)-(f), 
a = 0.6. 
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FIG. 2. Chaotic Spikes. The trajectory of U{0,t). Parame¬ 
ters: k = 3.7 X 10“^, /3 = 0.026999, and a = 0.6. 




FIG. 3. Lorenz maps of chaotic spikes, (a) The relation 
between and r^max- (b) The relation between t^^x and 

'^max- (c) The relation between t^ax and Simulation 

runs for 50000ts. Parameters are the same as in Fig. 


The above examples indicate that in the presence of 
STD, the dynamics of CANNs has a rich spectrum. The 
intrinsic dynamics is further enriched in response to ex¬ 
ternal stimuli. So, it is important to systematically study 
the structure of the dynamics of CANNs with STD. In 
this paper, we report our investigations of CANN with 
STD in two scenarios: one with very weak inhibition and 
the other with moderate inhibition in the presence of a 
single static input. 

For CANNs with STD and very weak inhibition, the 
dynamical picture is more complicated than that with 
moderate inhibition strengths, which only consists of 
phases with static bumps, metastatic bumps, moving 
bumps, and the silent state m- Depending on param¬ 
eters, the system can support the following behaviors: 
uniform firing patterns across the whole network, peri- 



FIG. 4. Phase diagram of CANN with STD. (a) Moderate 
global inhibition, (b) Weak global inhibition. Moving bump 
exists on the right side of the solid line. Light gray area 
(SAS) is for spikes and anti-spikes. Dark gray area (HS) is 
for homogeneous spikes. Dashed lines are boundaries of the 
simplihed dynamics (Eqs. and 0)- They are discussed 
in region B in Fig. The dot-dashed line is the boundary 
between the uniform hring state and the moving bump state 
given by Eq. |T^ . 


odic excitements of uniform firing, wave instabilities with 
complex and even chaotic firing patterns. 

For CANNs with STD and a single static input, the 
network can support various complex firing patterns. 
Those complex patterns are due to the interplay between 
the STD-driven intrinsic dynamics and the external in¬ 
put, depending on strengths of the single static input, in¬ 
hibition and STD. Similar to the previous scenario, some 
complex patterns are chaotic. 

The rest of the paper is organized as follows. We 
will begin with an introduction of the model being used 
throughout the paper. Next, we will present the intrinsic 
dynamics of the CANN with STD and very weak inhibi¬ 
tion. After that, we will report how the dynamics of the 
system is affected by the external input. Detailed behav¬ 
iors of the network in various phases in the phase diagram 
will also be reported. The relevance of the results will be 
discussed at the end. 


II. INTRINSIC BEHAVIORS UNDER WEAK 
GLOBAL INHIBITION 


N neurons are evenly distributed in the space of pre¬ 
ferred stimuli (in the following simulation results, N = 
256). Neurons are labeled by their preferred stimulus 
X. The range of {x} is (—1//2,1//2]. So the size of the 
space is L. Since usually the model is applied to the 
representation of directions or orientations, L = 27r and 
the periodic boundary condition is imposed. We mod¬ 
ify the general form of neural field theory and formulate 
the intrinsic dynamics of the neuronal input U{x^t) as 
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where p is the density of neurons in the space of preferred 
stimuli, I{x,t) is the external stimulus, J{x^x') is the 
coupling strength between neurons with preferred stimuli 
X and x', r(x', t) is the firing rate of neuron x' at time t. 
They are given by 
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( 3 ) 


where k is the strength of the global inhibition and 
[X]+ = max(X, 0). Here we adopt a Gaussian coupling 
and incorporate inhibitory connections into the global 
inhibition. 

p(x, t) is the availability of neurotransmitters in neuron 
X. The dynamics of p(x, t) is given by 


dpjx, t) 
dt 


Td 


/3p{x,t)r{x,t), 


( 4 ) 


where Td is the time scale of neurotransmitter recovery, 
which is chosen to be Td = SOrg. [3 is the fraction of to¬ 
tal neurotransmitters consumed by firing per spike. This 
implies that when the firing rate r is 0, p will gradu¬ 
ally recover to 1 and the effective connection strength 
will be exactly given by J{x,x'). However, p is less 
than 1 for active neurons and the connection strength 
is undermined by inefficient transmission. For simplic¬ 
ity of analysis, we introduce the rescaled variables and 
parameters: U = pJoU^ [3 = TdP/{pJoY^ ^ 
k = SV^ak^p^‘^ and I = pJo/. So that we could 
rewrite Eqs. ([H-0 as below. 
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In the absence of external inputs {I{x^t) = 0), a vari¬ 
ety of interesting behaviors have been discovered, such as 
the static bump, the moving bump and the silent state 
(Fig. [^a)-(c), see also USD- The static bump state, also 
known as a persistent spatially localized activity state 
[24] . is of interest because it can be found in physio¬ 
logical recordings in the prefrontal cortex during spatial 
working memory tasks and other systems that encode 
directional or spatial information, such as head direction 


cells in thalamus and basal ganglia and place cells in the 
hippocampus. The existence and stability of the static 
bump state were first analyzed in detail by Amari [2], 
followed by various extensions. Because of the transla¬ 
tional invariance of the neuronal coupling J(x,x'), the 
center of the static bump can be arbitrarily positioned. 
This is referred to as neutral stability, which leads to the 
naming of continuous attractor neural networks and the 
remarkable tracking ability of this model. 

The moving bump state corresponds to traveling waves 
which have been extensively studied experimentally ESl- 
im and theoretically [28II32] . In our model, neurotrans¬ 
mitters are depleted at the bump’s position due to the 
STD. Thus, the bump tends to move away to regions 
where neurotransmitters are more available, which is 
analogous to the spreading of forest fire as the fire also 
moves from places where trees are all burnt out, to places 
where fuels are abundant. This spontaneous movement 
is proposed to be able to compensate various kinds of 
delays in the neural system and therefore facilitate an 
accurate representation of moving stimuli m- However, 
these behaviors are found where the global inhibition k 
is relatively large. The parameter region where k is very 
small is not much explored yet. 


A. Phase Diagram 

The global inhibition plays an important role in shap¬ 
ing the bump states. Regions with low activity will 
be suppressed by regions with high activity through the 
global inhibition, thus the neuronal activity will be local¬ 
ized to form a state with only one bump, either moving 
or static. In the very weak global inhibition scenario, we 
expect that neural activity would be more uniform and 
synchronized [33]. In the uniform firing state, the firing 
rates of all neurons are uniform and time-independent 
(Fig-Bd)). 

Another possible effect of very weak global inhibition is 
that, neuronal activity will be relatively higher and more 
neurotransmitters are consumed, which induces popula¬ 
tion spikes. The neuronal activity grows up very fast 
due to the weak inhibition and a large amount of neu¬ 
rotransmitters are consumed. Then, the activity will die 
down because of the STD. This is called a population 
spike (also called ensemble synchronization) [34l|35]. Af¬ 
ter the activity dies down, neurotransmitters will grad¬ 
ually recover to the level that can support another pop¬ 
ulation spike. These population spikes show various dy¬ 
namics, including homogeneous spikes and the spike-and- 
anti-spike state, in the very weak global inhibition sce¬ 
nario (Fig. Be-f)). 

In the homogeneous spike state, all neurons in the net¬ 
work are synchronized in population spikes (Fig. Be)). 
In the spike-and-anti-spike state, one spike emerges at 
a certain place in the network, then splits into two sym¬ 
metric branches of moving bumps, which will collide with 
each other at the opposite side of the network and form 
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FIG. 5 . (Color online) Different behaviors of simplified ho¬ 
mogeneous dynamics. Variables are rescaled as u* = and 
p* = (p — 0 . 5 )^/ 0 . 5 ^ + 0.5 to make a readable presentation of 
trajectories while keeping topological equivalence. Gray area 
is the basin of attraction for the stable fixed point or the sta¬ 
ble limit cycle. Some trajectories are colored red to illustrate 
their approach to a stable limit cycle or stable fixed point, (a) 
Uniform firing without an unstable limit cycle. Parameters: 
/c = 5 X 10 “^, /3 = 0 . 025 . (b) Homogeneous spikes without an 
unstable limit cycle. Parameters: /c = 1 x 10 ~^, jS = 0 . 023 . 
(c) Homogeneous spikes with an unstable limit cycle. Param¬ 
eters: k = 2.9 X 10 “^, [3 = 0 . 0253 . (d) Uniform firing with an 
unstable limit cycle. Parameters: k — 7 x 10 -^ 13 = 0 . 0285 . 
(e) Silent state. Parameters: k = 3 x 10 “^, /3 = 0 . 027 . For 
(a)-(e), a = 0.6. 



FIG. 6. Phase diagram of simplified homogeneous dynamics 
for interaction range a = 0.6. Region A, B, C, D and E 
correspond to behaviors in Fig. [^a), (b), (c), (d) and (e), 
respectively. 


height of a spike is almost linearly dependent on the rest¬ 
ing time before that spike, rather than after that spike. 
Since longer resting times mean fuller recovery of neuro¬ 
transmitters, the population spike following a long rest¬ 
ing time is able to reach a greater height. 

Simulations are performed to find the boundaries of 
these behaviors in the phase space (Fig. [^. In the phase 
diagram, uniform firing phase is found where k is less 
than 4 X 10“^ and (3 is less than 0.02. When k increases, 
the uniform firing state becomes moving bumps. This 
is because stronger global inhibition will suppress the 
homogeneity of the dynamics. When p increases, the 
uniform firing state becomes homogeneous spikes (dark 
gray area in Fig. |^b)). The homogeneity of the dynam¬ 
ics is maintained, but stronger STD introduces temporal 
modulation. Spikes and anti-spikes (light gray area in 
Fig. I^b)) are found when P is even larger. Some part 
of this area overlaps with moving bump region, which 
indicates the coexistence of the two behaviors. 


B. Simplified homogeneous dynamics 


the anti-spike, and the firing activities stop. The cycle 
then repeats itself periodically (Fig. Bf))- 

In some parameter regions, spikes and anti-spikes can 
be chaotic. In the simulation, the initial condition is cho¬ 
sen to make spikes and anti-spikes appear at x = 0 and 
X = TT, respectively. We could see the chaotic behav¬ 
ior by just looking at the dynamics at x = 0 (Fig. [^. 
Define peak of t/(0,t) and as the 

temporal interval between the and (n + 1)^^ peaks 
of [7(0, t). By examining the relation between 

(Fig-Ba)), we could see it shows chaotic features. 
The relation between temporal intervals between spikes 
and the height of spikes (Fig.j^b) and (c)) shows that the 


In the phase diagram shown in Fig. [^b), there is a 
phase called uniform firing. To study the stability of this 
pattern, we simplified the dynamics in Eqs. and to 
be a pair of spatially independent differential equations: 


Td- 


du (t) 
^ dt 
dpjt) 
dt 


= -u (t) + 


p{t)u{ty 
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Ja^ 


= l-p{t)- 


Pp{t)u{t) 
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(8) 

( 9 ) 


Here, u{t) represents the spatially independent U{x,t). 
Ja = erf[L/(-\/8a)] and B = 1 + kL/{8'/^a)u^. The 
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nullclines for du/dt and dp/dt are 
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_ pu- 
^ ^ Ja-) 

Ppu^ 
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( 10 ) 


Intersections between nullclines are fixed point solutions. 
At the intersections, 


u = 


Ja A 


^ 8 \^a) 


( 11 ) 


solution with smaller the determinant of the matrix in 
Eq. ( p!^ i s negative. Hence it is a saddle point (triangles 
in Fig. [^. 

At the fixed point solution with larger the determi¬ 
nant is always larger than zero. The stability condition 
is given by the trace of the matrix. In order to have a 
stable fixed point solution, the trace should be less than 
zero. After some algebra, the phase boundary of uniform 
firing is given by 



L y'Td J (2-3)^ 

TsJI 2-{l + Ts/Td)B 

Td (2-Bf 


(13) 


Besides fixed point solutions, limit cycles could also 
appear through the Hopf bifurcation or the homoclinic 
bifurcation [36]. In Fig. there are typical phase por¬ 
traits in the simplified system for different values of k 
and p. In Fig. [^a), there are two stable fixed points, 
one representing the silent state and the other represent¬ 
ing the uniform firing state. The stable manifold of the 
saddle point is the boundary separating basins of attrac¬ 
tion of the two stable fixed points. Temporal oscilla¬ 
tions exist in the transient state, but the firing rate is 
time-independent at the steady state. The fixed point 
for uniform firing loses its stability in Fig. |^b), and a 
stable limit cycle appears around it, which demonstrates 
a Hopf bifurcation. The stable limit cycle corresponds 
to homogeneous spikes. Between Fig. |^a) and (d), a 
homoclinic bifurcation occurs. The unstable manifold of 
the saddle point touches its stable manifold, resulting in 
a homoclinic orbit, which is also an unstable limit cycle 
encircling the basin of attraction of the uniform firing 
fixed point. Beyond the homoclinic bifurcation, the un¬ 
stable limit cycle is no longer homoclinic and shrinks in 
Fig. [^d). Although this unstable limit cycle does not 
affect the stability of the fixed points, it confines the 
basin of attraction of the uniform firing state. Between 
Fig. I^b) and (c), the same homoclinic bifurcation hap¬ 
pens and the unstable limit cycle is the basin boundary 
of the stable limit cycle. Between Fig. [^c) and (d), the 
Hopf bifurcation same as the one between Fig. |^a) and 
(b) happens. Between Fig. |^c) and (e), the stable limit 
cycle and the unstable limit cycle approach each other 
and collide, which makes a fold bifurcation of cycles [36] . 

To study the stability of the above fixed point solu¬ 
tions, we assume that u{t) = uq ^ ui{t) and p{t) = 
Po + Pi {t ): where uq and po are the fixed point solution 
being investigated. By linearizing Eqs. (§ and (|^, we 
have 


where B > 1. This is a parametric expression for the 
boundary of the Hopf bifurcation (the dashed line is 
Fig.§. Other bifurcation curves in the phase diagram of 
this simplified system (Fig.[^ can be found using the nu¬ 
merical continuation package MATCONT m- We can 
see that this simplified model captures homogeneous be¬ 
haviors very well. However, for inhomogeneous behav¬ 
iors, we need to consider the wave stability of the fixed 
point solutions of Eqs. (§ and 


C. Wave stability of the uniform firing 


To understand why moving bump disappears in such 
low values of inhibition and what small parameter charac¬ 
terizes the existence of uniform firing, we consider fluctu¬ 
ations with wave vector g, so that U (x, t) = 
and p{x,t) = Po where uq and po are given 

by the steady state solution. Putting these into Eqs. 
and (§, and keeping terms up to the first order, we have 


dt 


Ui 

Pi 


2pouoQ-B 

TsB 

_ 2^poUq 
TdB 



(14) 


where Q = /tf /2 exp[-- iq{x - x')]- 

Note that here we just consider wave fluctuations, thus 
B is not affected by the fluctuations to the first order. 

To evaluate Q, we can express the integral in terms of 
Her mite polynomials. 




erf 
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iqa 

71 
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d f ui 
dt \ Pi 


2poUQJa—B Uo‘^ J a 

TsB TsB 

_ ‘2‘^PqUq _ B-\-^U(p‘ 

TdB TdB 



There are two non-zero fixed point solutions to Eqs. § 
and one with large u and the other with smaller as 
shown in Eq. 0 (see Appendix [a]). At the fixed point 


Owing to the smallness of the Gaussian factor in the 
higher order terms, we can approximate Q by Ja^~^ ^ ^^5 
where Ja = erf [l//(\/ 8 a)], or approximately 1 when a ^ 

L. 

For fixed point solutions, the stability condition is the 
trace of the stability matrix, T < 0. Hence the phase 
boundary is given by T = 2poUoJa^~^ ^ jB — \ — 

























6 


{l ^ Ts/Td = 0. Combining with the steady 

state solutions, we can express the stability condition as 

^ , (16) 

Numerical solutions show that the instability comes 
from the long wavelength mode. For a ring model of 
length L, the fundamental mode has a wave number 
q = 211 !L. Hence for L = 27r, g = 1. Combining with re¬ 
sults on the steady state solutions, the boundary of wave 
stability is 


k = 



Po 

1 -Po 


P- 


1 

(l-Pof 



(17) 


where the value of po on the boundary is — 

l)~^Ts/Td. This boundary is plotted as the dot-dashed 
line in Fig.j^b) and it agrees very well with the boundary 
between the uniform firing and the moving bump. The 
maximum of k on this boundary is 


k 


max 


2y/^J^ Ts‘^a 

(^2e-27r2aVL2 _ rd‘^L‘ 


(18) 


This shows that the smallness of k comes from the scaling 
{Ts/rd)^ a/L. On the other hand, the maximum of P on 
the boundary is 


/^max 


- 1 - T./Td jPs 


(2. 


,-27r2a2 


-ly 


Td 


(19) 


This shows that the smallness of (3 comes from the scaling 
TsiTd- 


III. RESPONSES TO A SINGLE STATIC INPUT 

In this section, we will discuss the network behavior in 
the presence of a single static input and moderate global 
inhibition. For the external input in Eq. we adopt 
the form /(x, t) = Aex.p[—{x — 2 ))^/( 2 a^)], where z is the 
center of that input (without loss of generality, ^ = 0 in 
this work), A is the strength of the input, and a a is the 
width of the input. Note that the behavior of this system 
is controlled by three parameters, namely the strength of 
the global inhibition k, the strength of the STD (3 and the 
strength of the external input A. In this work, different 
response patterns are discussed in the parameter space 
spanned by these three quantities. 



FIG. 7. (Color online) Four basic dynamic responses to a 
single static input in CANN with STD. The color scale shows 
the firing rate r{x,t). (a) Emitter. Parameters: k = 0.2, 

[3 = 0.3. (b) Population spikes. Parameters: k = 0.3, [3 = 0.4. 
(c) Moving bump. Parameters: k = 0.3, /3 = 0.1. (d) Slosher. 
Parameters: k = 0.5, [3 = 0.1. For all (a)-(d), A = 0.8, 
a = UA = 48° = 0.8378. 


(a) 




Mixture behavior 
Bistable 




k 
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FIG. 8. Phase diagrams for the four basic responses in the 
space of k and [3 with different values of A. a = 0.5 and 


aA = \/2/2. 


A. Four Basic Dynamic Response Patterns 

The static bump is expected to be the simplest form of 
response. However, in a very large region of the parame¬ 
ter space, static bumps are unstable and much more in¬ 
teresting response patterns emerge. Among them, there 


are four basic dynamic patterns through which we can 
understand the general property of this system (Fig.I^a)- 

(d)). 

One response pattern is the moving bump. Moving 
bumps result from the mobility of the neural field en¬ 
hanced by the STD. Once a bump is built, neurotrans- 
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FIG. 9. (Color online) Phase diagram in the space of A 
and (3. The color code indicates lengths of the periods in log 
scale. Black dots means the period of the response is too 
long to be detected by the program, suggesting the response 
is possibly chaotic. The global inhibition strength k — 0.3. 
a — aA — 0.8378. (a) Full model simulation. M, E, P and 
S represent “moving bump”, “emitter”, “population spikes” 
and “slosher”, respectively, (b) Numerical solutions of the 
second order Fourier series expansion. The gray regions are 
bistable regions where the lengths of the periods can be ei¬ 
ther of two different values. The box in (b) encircles the 
region where Lyapunov exponents are computed and shown 
in Fig. 
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FIG. 10. (Color online) Largest Lyapunov exponents (LLE) 
computed on the equations of second order Eourier series ex¬ 
pansion. Parameters are the same as those in Eig. E Gray 
scale denotes the lengths of the periods for periodic behaviors. 


mitters are depleted in the bump region, leading to a 
tendency of the bump to move away to fresher regions. 
As an intrinsic behavior, the moving bump will keep its 
profile and its speed all the time (Fig. [^c)). However, 
while the static input is imposed, the speed and profile 
of the bump will change when the bump crosses the in¬ 
put. The bump is higher and faster when approaching 
to the input, while weaker and slower when leaving the 
input, because the external input tends to attract the 
bump (Fig.j^c)). 

When the attraction provided by the external input 
is strong and the mobility enhanced by the STD is not 
sufficient for the bump to overcome the attraction of 
the input, the bump gets trapped and moves side-to- 
side around the external input. It is called a slosher 

(Fig. gel)) m- 

In both cases of the moving bump and the slosher, the 
dynamics are governed by the mobility enhanced by the 
STD and the attraction provided by the external input. 
The amplitude of the bump does not change significantly 



FIG. 11. (Color online) Phase diagrams in the space of A 
and 13 with different values of k. k is 0.25 in (a), 0.35 in (b), 
0.40 in (c), 0.45 in (d), 0.60 in (e) and 0.70 in (f). The color 
code indicates the lengths of the periods in log scale. Black 
dots mean that the period of the response is too long to be 
detected by the program, suggesting the response is possibly 
chaotic. M, E, P and S represent “moving bump”, “emitter”, 
“population spikes” and “slosher”, respectively, a = aa = 
0.8378. 



EIG. 12. (Color online) Phase diagram in the space of k and 
[3. The color code indicates the lengths of the periods in log 
scale. Black dots means that the period of the response is too 
long to be detected by the program. Parameters: A = 0.8, 
a = 0.5 and ua = V^l2. 


during its movement. However, when p and A are suf¬ 
ficiently large, the effect of the STD is not just mobility 
enhancement, but also amplitude modulation. Large [3 
and A means that the bump will consume more neuro¬ 
transmitters so that it cannot maintain its amplitude all 
the time. The emitter (Fig. [^a)) is an example in such 
case. One moving bump is emitted by the external input. 
After it travels around the network, the bump dies down 
due to the excessive consumption of neurotransmitters 
during traveling. Then, the network waits a while un¬ 
til sufficient amount of neurotransmitters is recovered to 
support another emission of the moving bump. When 
































FIG. 13. (Color online) Mixtures of emitters and population 
spikes, (a) Periods of different responses with different input 
strength A. Dots are simulation results. If the period is too 
long to be detected, there is no dot. Six open circles are places 
where examples in (b)-(g) are drawn from, (b)-(g) Examples 
of mixture behaviors between emitters and population spikes. 
They correspond to the circles in (a) from left to right. Color 
scale indicates the firing rate r(x,t). Each of (b)-(g) shows 
one period of the particular behavior, k = 0.3, /3 = 0.3, and 
a = aA = 0.8378. 


EIG. 14. (Color online) Mixtures of moving bumps and 
sloshers. (a) Periods of different responses with different input 
strength A. Dots are simulation results. If the period is too 
long to be detected, there is no dot. Six open circles are places 
where examples in (b)-(g) are drawn from, (b)-(g) Examples 
of mixture behaviors between moving bumps and sloshers. 
They correspond to the circles in (a) from left to right. Color 
scale indicates the firing rate r(x,t). Each of (b)-(g) shows 
one period of the particular behavior, k = 0.3, 13 = 0.13, and 
a = aA = 0.8378. 


the external input is even stronger, we see a similar re¬ 
sponse, namely, population spikes (Fig. [^b)), in which 
case a static bump, rather than a branch of moving bump, 
is emitted after recovery, since the external input is so 
strong that the bump is trapped [34l [35] . This behavior 
is similar to breathers [HI man, except that breathers 
oscillate in their widths, whereas population spikes pri¬ 
marily oscillate in their heights. 

We explore the four basic dynamic responses in the 
parameter space of k and P with different values of A 
(Fig. §. When the input is weak (Fig. [^a) and (b)), 
sloshers appear between the moving bump region and 
the static bump region, compared to the intrinsic behav¬ 
iors in Fig. [^a). This is because the mobility enhanced 
by the STD is not enough to delocalize the bump, which 
is attracted by the static input. When the input is strong 
(Fig. I^c) and (d)), the consumption of neurotransmit¬ 
ters is so fast that the bump cannot keep its amplitude 
stable. This results in the emergence of the emitter and 
population spikes. We also notice that when the input 
is weak, there is a bistable region for static bumps. In 
this bistable region, there are two stable static bump so¬ 
lutions, corresponding to the self-sustained bump and a 
weaker bump that is created only when M > 0. A bifur¬ 
cation diagram of these static bump solutions is shown 


in Fig. jT and discussed in Appendix [B] 


B. Mixture Behaviors 

In numerical solutions, we find that in a large part of 
the parameter space, response patterns can be none of 
the four basic patterns and very complex. They seem 
to be different mixtures of the four basic dynamic pat¬ 
terns. Most of these responses are periodic. The tempo¬ 
ral duration of one period is closely related to what kind 
of responses are being mixed. Similar behaviors have 
been observed in other models of neural fields with spike 
frequency adaptation or short-term synaptic depression 
[40[|4T]. However, relations between different mixture be¬ 
haviors have not been systematically understood. Here, 
we proposed that this can be done by monitoring the 
period of asymptotic states. Thus, we can show how 
different mixture behaviors are organized in the phase 
diagram in the space of A and [3 (Figs. and [TT|, and in 
the space of k and P (Fig. [T^. 

For the full model simulation (Fig. [^a)), where the 
number of neurons is 256, Eqs. 0 are solved by us¬ 
ing the MATLAB command ode45 and the period is de¬ 
termined by examining the auto-correlation function of 
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FIG. 15. (Color online) Mixtures of emitters and population 
spikes, (a) Periods of different responses with different global 
inhibition strength k. The vertical axis is in logarithmic scale. 
Dots are simulation results. If the period is too long to be de¬ 
tected, there is no dot. Nine open circles are places where 
examples in (b)-(j) are drawn from, (b)-(j) Examples of mix¬ 
ture behaviors between emitters and population spikes. They 
correspond to the circles in (a) from left to right. Color scale 
indicates the firing rate r(x,t). Each of (b)-(j) shows one pe¬ 
riod of the particular behavior. /3 = 0.36, A = 0.8, a = 0.5 
and a A = a / 2 / 2 . 


EIG. 16. (Color online) Mixtures of moving bumps and slosh- 
ers. (a) Periods of different responses with different global in¬ 
hibition strength k. The vertical axis is in logarithmic scale. 
Dots are simulation results. If the period is too long to be 
detected, there is no dot. Nine open circles are places where 
examples in (b)-(j) are drawn from, (b)-(j) Examples of mix¬ 
ture behaviors between moving bumps and sloshers. They 
correspond to the circles in (a) from left to right. Color scale 
indicates the firing rate r(x,t)- Each of (b)-(j) shows one pe¬ 
riod of the particular behavior. /3 = 0.1, A = 0.8, a = 0.5 and 
aA = ^/2/2. 


asymptotic states. In Fig.j^b), the period is determined 
by solving a boundary value proble m o f the second order 
Fourier series expansion (see Eq. (B3) in Appendix H 
while letting A increases or decreases, using the MAT- 
LAB command bvp4c. 

In the phase diagram, there are many patches within 
which the period of the dynamics changes continuously. 
However, the period jumps abruptly across boundaries 
of the patches. This indicates that behaviors are sim¬ 
ilar within each patch, while transitions happen across 
boundaries. Boundaries are relatively coarse in Fig. ©a) 
due to numerical errors and the bi-stability of the dynam¬ 
ics, which is shown clearly in Fig.j^b). In the gray region 
along phase boundaries, different dynamics can be found 
by starting with initial conditions from different sides of 
the boundaries. The four basic dynamic response pat¬ 
terns are located at the four disjoint regions of the phase 
diagram. In between them, there is a rich spectrum of 
different mixture behaviors. Especially, black dots, where 
the length of the period is too long to be well determined 
within a time limit, are found in the mixture behavior 
region, which may imply chaos. Largest Lyapunov expo¬ 
nents are computed using Wolf’s algorithm [42l|43] in re¬ 
gions containing black dots (Fig. 10). Positive exponents 


exist extensively around /3 = 0.2, which clearly demon¬ 
strates the existence of chaos in this strongly coupled 
neural field. In the green area to the left of the sloshers 
region, sloshers are quasi-periodic due to the Neimark- 
Sacker bifurcation [36] illustrated in Fig. 17 in Appendix 
[B| The similarity between Fig. |^a) and (1^ justifies the 
method of using the Fourier basis to study this system. 


C. Phase Diagrams 

Having explored the phase diagrams in different situa¬ 
tions, we find that the parameter space can be separated 
into two parts. The upper part where /3 is larger consists 
of emitters, population spikes and their mixtures (Fig.[T^ 
and 


15). The lower part where [3 is smaller consists of 


moving bumps, sloshers and their mixtures (Fig. [Mjand 
16). Between these two parts around p = 0.2, responses 
tend to have very long period and show chaotic features. 
We may conclude that short-term synaptic depression 
have different effects depending on its strength. 

Weak STD enhances the mobility of the bump without 
affecting the amplitude of the bump significantly, leading 
to bumps of relatively stable amplitude and varying posi- 
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tion. This is the spatial modulation effect of the STD. In 
Fig. we see that as the external input gets stronger 
and stronger, the mixture behaviors tend to have less 
and less emitter components but more and more slosher 
components, from Fig, p^b) to (g). 

On the other hand, strong STD disrupts the bump in 
time, since neurotransmitters are depleted rapidly during 
the spikes. This shows the temporal modulation effect of 
the STD. Bumps in the time sequence generally are not 
the same, implying a possibility to encode different in¬ 
formation in different emissions, an example having been 
discussed in detail in m- In Fig.f^ we see that emission 
of population spikes progressively dominate over emission 
of moving bumps as A increases from Fig, p^b) to (g). 

Varying the global inhibition at a constant external 
stimulus strength shows similar effects. Stronger inhibi¬ 
tion suppresses activities in the region outside the exter¬ 
nal stimulus, therefore confines the bump or spikes to the 
region near the external stimulus. Hence for the case of 
strong STD shown in Fig.[^ the emitter components in 
the mixture behavior are progressively replaced by popu¬ 
lation spikes when global inhibition increases. Similarly, 
for the case of weak STD in Fig. the moving bump 
components in the mixture are progressively replaced by 
sloshers when global inhibition increases. 


IV. DISCUSSION 

We have found a rich spectrum of firing patterns in 
CANNs with STD in the regime of weak inhibition and 
the regime of a single static input. CANNs with moder¬ 
ately strong inhibition were initially introduced to track 
continuous inputs using static and moving bumps. How¬ 
ever, in the very weak global inhibition region, CANNs 
can no longer support bumps. Instead, the dynamics of 
the network show population spikes of various kinds. In 
particular, chaotic behavior is found in the amplitudes 
of the population spikes, typical of cycles of storing and 
releasing resources (neurotransmitters) in pulses. The 
smallness of STD scales as and the smallness of the 
global inhibition scales as {ts/T d)^{a/L). 

In the case of a single static input, we have found four 
basic patterns of dynamic responses and their mixtures. 
When STD is weak (lower than 0.2, roughly), it mainly 
provides spatial modulation, or in other words, enhances 
the mobility of the bump. Inputs of different strengths 
provide different attraction, leading to moving bumps, 
sloshers or mixtures of them. When STD is strong, STD 
provides temporal modulation along with spatial modu¬ 
lation. Together with the static external input, they re¬ 
sults in emitters, population spikes, or mixtures of them. 
In the parameter region where STD strength is interme¬ 
diate, chaotic behaviors appear. Although it is not fully 
understood, we believe that the involvement of both tem¬ 
poral and spatial modulation of STD is the major cause 
of complexity in that region. 

Due to their richness, the firing patterns have poten¬ 


tials in encoding and decoding information. An example 
is the decoding of two inputs that fluctuate in time and 
overlap in the space of preferred stimuli so strongly that 
their time average becomes indistinguishable. Hence any 
time-independent decoding methods are rendered ineffec¬ 
tive. In m we demonstrated that temporal modulations 
of the population spikes can provide a mechanism to re¬ 
solve the two inputs, and produce results consistent with 
the resolution enhancement in transparent motion m- 

Bifurcation analysis has provided important insights 
into the neural field dynamics, indicating underlying 
mechanisms for a variety of dynamical behaviors. In mi, 
a similar neural field model with linear spike frequency 
adaptation was studied with detailed bifurcation analy¬ 
sis. In the presence of a simple weak input, their model 
also shows static bumps, moving bumps, and sloshers. 
The population spikes, which need a relatively strong in¬ 
put (Fig. 1^, were not found in [44] because the inputs 
were relatively weak. A potential mechanism for percep¬ 
tion switching with complex inputs mediated by sloshers 
was proposed. In [45], a mechanism leading to chaos was 
shown in effect in the single neuron model with short¬ 
term synaptic plasticity. The chaotic behavior emerges 
through a Shil’nikov bifurcation of homoclinic orbits. To¬ 
gether with our mechanism on the population level, how 
the short-term synaptic plasticity enhances signal pro¬ 
cessing in the neural system across multiple levels should 
be of interest for further studies. 

Furthermore, the existence of chaotic behaviors in our 
network is relevant to the so-called “edge-of-chaos” re¬ 
gion, which has been observed to coincide with best com¬ 
putational competence Hi HZ]. Near the edge of chaos, 
the behavior of the network is neither dominated by the 
internal dynamics so as to be insensitive to external in¬ 
puts, nor does it depend on the external perturbations 
so much as to be vulnerable to any noise. Although 
some argued that operating near the edge of chaos is nei¬ 
ther a sufficient nor necessary condition for the system to 
achieve the best computational power [48] , chaotic neural 
networks are still often resorted to in modeling generic 
cortical microcircuits [49] . 

Our work shows that even a recurrent network with 
a highly regular structure can support extremely com¬ 
plex dynamics and chaos, in the presence of short-term 
synaptic plasticity. Previous work showed that when 
randomly connected recurrent neural networks exhibit 
chaotic activities, they act as a dynamical repertoire pow¬ 
erful in performing a variety of complex computational 
tasks [50||5T] and capable of reproducing main features of 
certain experiments [52] through different learning rules. 
However, the randomness in their connectivity and the 
overwhelming richness of their dynamics make theoreti¬ 
cal understanding and predictable generalization difficult 
[531154], although in some cases dynamical skeletons can 
be extracted to elucidate how the network achieves dif¬ 
ferent functions [55]. In CANNs with STD, the edge of 
chaos emerges in the region where the primary effect of 
STD changes from spatial modulation to temporal mod- 
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ulation. A rich spectrum of dynamical behaviors can 
be readily found and understood near the edge of chaos. 
How to tap into the potential computational power of 
CANNs with STD is an interesting problem to be inves¬ 
tigated in the future. 
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Appendix A: stability analysis of fixed points in the 
simplified homogeneous dynamics 


The trace T and determinant D of the Jacobian matrix 
of Eqs. ^ and are, 

/3uo‘^' 


Td 


B 


D = 


Ts{U{)Ja 2 ) 
TdB 


For the two fixed point solutions in Eq. (11), 


D = ±- 


(Al) 

(A 2 ) 

(A3) 


where. 


1 

ui{t) = — u{x^t)e 

L J_LI2 

1 

{t) = - (B2) 

B J-L/2 


Pi{t)=L 

Therefore, the CANN model can be rewritten as 


0,2 7r2/2/r2 




^ ^ a-s^v=l 


q—s-\-v=l 


^ ' q-s+v=l 


B(t) — 1 H- 


kL 


8\/^a 


M 

E (B3) 

= -M 


where Si^o is equal to 1 only when / = 0 , otherwise 0 . 
14*(t) is the complex conjugate of Us{t). g, s, v are all 
integer indices ranging from —M to M. 

Eor M = 0, Eq. (B3) reduces to the simplified ho¬ 
mogenous dynamics in Eqs. ^ a nd Eor M = 1, 
Eq. (B3) is equivalent to Eq. (14) where the wave sta¬ 


bility is analyzed. Eor M = 2, it would be complicated 
to apply analytical methods. However, numerical meth¬ 
ods are far more efficient for the Eourier series expansion 
up to the second order than for the full network model. 
Moreover, almost all features of the phase diagram of the 
full net wor k model are maintained in the phase diagram 
of Eq. (B3) with M = 2 (see Eig. [^a) and (b)). 


where 7 = /3 + kL/{sV^a) and plus signs are for the 
solution with larger u. D is the product of the two eigen¬ 
values. Therefore, the fixed point with smaller 14 is a 
saddle point with D < 0 (triangles in Eig. [^, whereas 
the stability of the fixed point with larger u [D > 0) de¬ 
pends on the sign of the trace T. Combining Eq. (Al) 


with the expression for the nullclines in Eq. ( |1Q| ) and the 
definition for H, we can derive the parametric expression 
(Eq. ( p!^ ) for the boundary where T = 0. 


(a) 



A 


(b) 
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Appendix B: Fourier series expansion of the CANNs 

model 

For the CANN model Eqs. ([^-0 and the external 
input /(x,t) = Aexp[—(x — z)^/{2a^)], where a a is the 
width of the input, we expand 1/(x, t) and p(x, t) in terms 
of Eourier series up to order, 

M 

U{x,t)= E 

l=-M 

M 

p(rr,t)= ^ (Bl) 

l=-M 


FIG. 17. (Color online) (a) The real part of the first order 
component ui of solutions of Eq. (B3) with different strengths 
of external inputs A. The solid line is fixed point solutions. 
The red dot and blue dot labeled LP (limit points) corre¬ 
sponds to two saddle-node bifurcations, respectively. The 
black dot labeled H indicates a Hopf bifurcation. The dashed 
lines are extrema of limit cycles. The black dot labeled NS 
indicates a Neimark-Sacker bifurcation. Other parameters: 
k = 0.55, /3 = 0.08, a = 0.5 and ga = V^l2. (b) Continu¬ 
ation of the saddle-node bifurcation in the parameter space. 
The left and right surfaces corresponds to the red and blue 
dots in panel (a), respectively. The color scale shows the real 
part of the first order component ui. The cusp bifurcation 
happens along the intersection of the two surfaces. 
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A bifurcation diagram is computed numerically using 
MATCONT [37] to illustrate the bistable region in Fig.|^ 
and the emergence of sloshers. Within the region in 
Fig.jl^b), there are three static bump solutions centered 
at X = 0. The lower one and higher one corresponds to 
the weak bump that is created only when A > 0 and 
the self-sustained bump. The middle one is a saddle 


point which separates basins of attraction of the other 
two bump solutions. The slosher appears after a Hopf 
bifurcation of the higher static bump state. When A de¬ 
creases, it will become unstable after a Neimark-Sacker 
bifurcation [36| , which produces quasi-periodic behaviors 
shown in the green area to the left of the sloshers region 
in Fig. [TOj 
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